Gravitational Waves from the Dynamical Bar Instability in a Rapidly Rotating Star 
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A rapidly rotating, axisymmetric star can be dynamically unstable to an m = 2 "bar" mode that 
transforms the star from a disk shape to an elongated bar. The fate of such a bar-shaped star is 
uncertain. Some previous numerical studies indicate that the bar is short lived, lasting for only 
a few bar-rotation periods, while other studies suggest that the bar is relatively long lived. This 
paper contains the results of a numerical simulation of a rapidly rotating 7 = 5/3 fluid star. The 
simulation shows that the bar shape is long lived: once the bar is established, the star retains this 
shape for more than 10 bar-rotation periods, through the end of the simulation. The results are 
consistent with the conjecture that a star will retain its bar shape indefinitely on a dynamical time 
scale, as long as its rotation rate exceeds the threshold for secular bar instability. The results are 
described in terms of a low density neutron star, but can be scaled to represent, for example, a 
burned-out stellar core that is prevented from complete collapse by centrifugal forces. Estimates for 
the gravitational-wave signal indicate that a dynamically unstable neutron star in our galaxy can 
be detected easily by the first generation of ground based gravitational-wave detectors. The signal 
for an unstable neutron star in the Virgo cluster might be seen by the planned advanced detectors. 
The Newtonian/quadrupole approximation is used throughout this work. 
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I. INTRODUCTION 

A self-gravitating, axisymmetric fluid body with a suf- 
ficiently high rotation rate can be dynamically unsta- 
ble to non-axisymmetric perturbations. Typically, the 
fastest growing unstable mode is the m = 2 "bar" mode 
which acts to transform the body from a disk-like shape 
to an elongated bar that tumbles end-over-end. This in- 
stability has been described analytically for the case of 
uniform density bodies [|l],D, and has been the subject 
of numerous numerical studies^. The numerical results 
show that bar formation is accompanied by the ejection 
of mass and angular momentum, and that the ejected 
matter forms long spiral arms in the equatorial plane. 
The subsequent evolution is less certain. Some simula- 
tions indicate that the bar shape is short lived, with the 
star returning to a predominantly disk-like shape after a 
few bar-rotation periods. Other simulations predict that 
the bar persists for many bar-rotation periods. In recent 
work, New, Centrella, and Tohline address this issue 
with a series of simulations using two different codes at 
various resolutions, and conclude that the bar shape is 
persistent. In their highest resolution run the bar de- 
cayed after roughly 6 or 7 bar-rotation periods. This 
was believed to be caused by numerical errors that in- 
duced an unphysical center of mass motion. In a lower 
resolution run a symmetry condition was imposed that 



prevented any center of mass motion. In that case the 
star maintained its bar shape throughout the simulation. 

The purpose of the present work is to simulate the 
long-time evolution of a rapidly rotating, self-gravitating 
star using Newtonian hydrodynamics and gravity. As 
discussed in Sec. [II, the numerical code is substantially 
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different from the codes that have been used previously 
to address this problem. The initial data for this study 
consists of a 7 = 5/3 polytrope with stability parameter 
fi = 0.30. The stability parameter is defined by /3 = 
T/|VF|, where T is rotational kinetic energy and W is 
gravitational potential energy. The results here suggest 
that the bar shape is indeed long-lived — the star displays 
a prominent bar shape at the end of the simulation, which 
includes more than 10 bar rotation periods. 

Numerical studies of fluids with various equations of 
state and initial rotation profiles have shown that the 
dynamical bar instability appears when the stability pa- 
rameter j3 exceeds a certain critical value, typically close 
to 0.27 For the secular instability, which arises 

through dissipative mechanisms, the critical value of (3 is 
near 0.14. A neutron star might reach the critical value of 
(3 for dynamical or secular instability by accreting mat- 
ter and angular momentum from a binary companion. 
A stellar core that has exhausted its nuclear fuel might 
reach a critical rotation rate as it collapses. 

A star or stellar core that develops a rotating bar-like 
configuration will generate large amounts of gravitational 
radiation. Depending on the distance of the source, this 
radiation might be strong enough to be detected by the 
world-wide network of gravitational-wave detectors cur- 
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rently under construction g]. Here, the question of per- 
sistence of the bar shape becomes very important. The 
detectabihty of a source depends on its characteristic am- 
pHtude he ~ hy/n, where h is the amphtude of the waves 
with frequency / and n is the number of wave cycles in a 
bandwidth near / Thus, long-duration signals with 
large n can be more easily detected than short-duration 
signals. 

Should we expect the bar shape to persist or decay? 
One reason why the bar might decay is the loss of mass 
and angular momentum from the ends of the bar. The 
accompanying drop in rotational kinetic energy could re- 
duce the stability parameter and allow the star to re- 
turn to axisymmetry. Loss of rotational kinetic energy 
through shock heating might also occur. For the sim- 
ulation presented here (3 has an initial value of 0.30, 
large enough to dynamically trigger the growth of the 
bar mode. During the initial period of bar formation, 
mass and angular momentum are shed from the ends of 
the bar and long, spiral arms are formed. The stabil- 
ity parameter rapidly drops below 0.27 and eventually 
settles to a value of about 0.24. However, these losses 
of mass and angular momentum, with the accompanying 
drop in /3, are insufficient to completely rob the star of 
its bar shape. 

The results here are consistent with the conjecture that 
the star will retain its bar shape indefinitely on a dynam- 
ical time scale, as long as /3 stays above the critical value 
for secular instability If this is correct, the bar should 
still decay due to secular processes. Imamura, Durisen, 
and Pickett have argued that the bar continues to shed 
small amounts of angular momentum to the spiral arms, 
and these loses cause the bar to decay |8j. Other secular 
mechanisms, such as energy and angular momentum loss 
through gravitational radiation, could play a role as well. 

It is possible that, for the equation of state and angular 
velocity profile used here, the star is able to retain its bar 
shape simply because the critical value of /3 for dynamical 
instability is 0.24 or less. This possibility was tested by 
conducting a second simulation using the same equation 
of state and scaling the angular velocity by a factor of 
3/4. The resulting model has an initial value of (3 equal 
to 0.25. The bar mode showed no signs of growth in this 
simulation. Thus, for the models tested here, the bar 
mode does not spontaneously grow unless (3 exceeds a 
critical value greater than 0.25. However, once the bar 
shape is established, it can persist for many bar-rotation 
periods with (3 equal to 0.24 or less. 

In the present numerical calculation there are two 
sources of error in angular momentum. First, some of 
the angular momentum that is lost off the edge of the 
computational grid might, ideally, fall back onto the star. 
Second, numerical errors at each time step introduce a 
purely artificial loss of angular momentum. Both of these 
errors act to decrease the angular momentum of the star. 
One expects such losses to cause the bar shape to de- 
teriorate. What is seen, both in the gravitational wave 
analysis and in the Fourier analysis of the density, is a 



relatively persistent signal with a modest decline in am- 
plitude over the duration of the simulation. 

The physical model and initial data that form the basis 
of the simulation are described in Sec. II. Section HI con- 
tains a discussion of the numerical code. The results of 
the simulation are described qualitatively and quantita- 
tively in Sec. IV. Conclusions and results are summarized 
in Sec. V. Details of the initial data code are given in the 
Appendix. 

II. PHYSICAL MODEL 

The initial data for the simulation consists of a station- 
ary, axisymmetric Newtonian fluid star with polytropic 
equation of state, P = Kp"'. The equation of state pa- 
rameters are chosen to coincide roughly with those of 
a low density neutron star with soft equation of state, 
7 = 5/3 and K = 5.38 x lO^cgs [§. The bar instabil- 
ity in stellar cores with stiff equations of state has been 
investigated by Houser and Centrella 

The equations of hydrostatic equilibrium are solved us- 
ing an algorithm described in the Appendix, in which the 
freely specifiable data are the central density pc and the 
angular velocity distribution w(r). Here, r is the distance 
from the rotation axis. With the choices 

Pc = 2.00 X lO^'^g/cm^ , (la) 
Loir) = (4000/s)e-(''/4-80xio'--)' , (lb) 

the equilibrium configuration has mass, equatorial ra- 
dius, polar radius, angular momentum, and stability pa- 
rameter given by 

M = 2.37 Mq, (2a) 

i?cq = 4.91 X 10** cm , (2b) 

i?p = 1.11 X lO^cm , (2c) 

J = 6.98 X lO^'^gcmVs , (2d) 

(3 = 0.300 , (2e) 

respectively. Equations (lb) and (2b) show that the an- 
gular velocity has values 

Lu{0) = 4000/s , t^(i?oq) = 1400/s , (3) 

on the rotation axis and at the equator. The azimuthal 
velocity at the equator is 6.89 x 10® cm/s, below the Ke- 
pler velocity of 8.90 x 10® cm/s. Note that the mass (2a) 
is considerably larger than the masses of observed neu- 
tron stars, and exceeds the limit for nonrotating neutron 
stars for most equations of state. However, as discussed 
by Baumgarte, Shapiro, and Shibata [|lO|, differentially 
rotating "hypermassive" neutron stars with M > 2 Mq 
might appear as the immediate products of core-collapse 
supernovae or binary neutron star mergers. 

The rotation law for a uniform density Maclaurin 
spheroid is 
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jsp(men) = (5J/2M)[1 - (1 - men/Af)2/3] 



(4) 



where jsp = r ui is the specific angular momentum 
and men is the mass enclosed in a cylindrical radius. 
Many previous works on the dynamical bar instability 
have em ploy ed this same rotation law (see, for example, 
Refs. |Tl ,p^,|9|,p| ) . For the initial data code used here the 
angular velocity is specified as a function of radius r, so 
the rotation law (4) cannot be used directly. However, 
the rotation law (lb) was chosen for its similarity to the 
Maclaurin law, as Fig. 1 shows. In Fig. 1, the munerical 
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FIG. 1. The specific angular momentum jsp, in units of 
J/M, versus the mass men enclosed in a cylindrical radius, 
in units of Af. Both the numerical data (crosses x) and the 
Maclaurin rotation law (solid curve) are shown. 



For the model with the scaling displayed in Eqs. (1) 
and (2), the maximum flow velocity is 8.2 x lO^cm/s. 
This occurs at a distance of 3.4 x 10^ cm from the rota- 
tion axis. At this peak velocity, the special relativistic 
gamma factor (1 — w^/c^)"^/^, where c is the speed of 
light, differs from unity by about 4%. General relativis- 
tic effects should be greatest near the surface of the star 
at the poles, where the gravitational potential is approx- 
imately GM/Rp = 2.8 x lO^^cm^/s^. Comparing this 
result with , we find that general relativistic effects rep- 
resent a correction of roughly 30% at the poles. Thus, 
with the scaling used in Eqs. (1) and (2), the errors that 
we introduce by ignoring special and general relativity 
are substantial. 

The errors are reduced if the model is rescaled ap- 
propriately. For example, with A = 1/2, /x = 2, and 
T = 1/4, the mass of the star is cut in half, its linear 
size is increased by a factor of 2, and its angular veloc- 
ity is decreased by a factor of 4. In this case the model 
might represent a stellar core that has partially collapsed, 
but is prevented from contracting to nuclear density by 
centrifugal forces. (Note that with this scaling New- 
ton's constant is unchanged. Also note that the speed 
of light should not be scaled, since it does not appear in 
the Newtonian hydrodynamical equations or the Poisson 
equation.) The scaling reduces the maximum speed by 
a factor of 2 and reduces the gravitational potential by 
a factor of 4. For this model, special and general rela- 
tivistic effects represent corrections of about 2% and 8%, 
respectively. 



data for jsp and men are plotted for the rotation law (lb), 
along with the Maclaurin rotation law (4). For clarity, 
only one in five numerical data points is shown. The 
data and the rotation law (4) agree on the value of jsp 
to within 2% throughout most of the star, out to about 
nien/M = 0.8. At the surface of the star, the rotation 
law (lb) differs from the Maclaurin rotation law by about 
8%. 

The dimensionful model described here can be rescaled 
by introducing parameters A, /x, and r that scale the 
lengths, masses, and times. To be precise, a quan- 
tity Q with dimensions [Q] = L^M^T^ , where L is 
length, M is mass, and T is time, is rescaled accord- 
ing to Qncw — Qoid/i^^^J-^T'^)- A physical, dimensionful 
model is retained if the parameters are dimensionless and 
Newton's constant G is unchanged by the rescaling. This 
requires A, /i, and r to satisfy A^ = ^t^. Alternatively, 
the model can be converted to polytropic units [n2| in 
which G = K = M ^ 1 hy setting A = 4.81 x lO^cni, 
fi = 4.71 X 10^3 and r = 1.88 x 10"^ s. This leads to 

Pc = 4.73 X 10"^ , (5a) 
i?cq = 10.2 , (5b) 
w(0) = 7.53 X 10"^ , (5c) 

for the central density, equatorial radius, and central an- 
gular velocity. 



III. NUMERICAL CODE 

The axisymmetric initial data was calculated on a grid 
in the r-z plane consisting of 512 x 511 zones. Details 
and tests of the initial data code are contained in the 
Appendix. In preparation for the evolution of this data, 
the grid was chosen to cover a physical domain of 1.77 x 
10^ cm in radius and 2.50 x 10^ cm in the z-direction. 
Thus, the star was contained in a small region of the 
grid, roughly the inner 142 x 45 zones. The quality of the 
data can be checked with the diagnostic expression V = 
{2T + W + 3J PdV)/W, where T is kinetic energy, W is 
gravitational potential energy, and J P dV is the volume 
integral of pressure. Accordin g t o the virial equations, 
this quantity should vanish ||l|,Q. For the numerically 
generated model, this diagnostic has a value of V = 3.2 x 
10"^ 

At the beginning of the evolution, the density and ve- 
locity arc interpolated onto a three-dimensional Carte- 
sian grid. The grid contains 128"^ zones and covers a cubic 
domain with sides of length 2.50 x 10^ cm. The equato- 
rial radius of the star spans 25 zones, while the polar 
radius spans 6 zones. The density in each zone is modi- 
fied with a random perturbation ranging from —10% to 
+10% of the unperturbed value. Initially, the specific 
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internal energy e is obtained from the density by the re- 
lation e = Kp^^^ ji^^— 1). Thereafter, the star is evolved 
with the equation of state 



P=(n-\)pe 



(6) 



with 7 = 5/3. 

The evolution code includes the nonrelativistic hy- 
drodynamics code VH-1, written by Blondin, Hawley, 
Lindahl and Lufkin. VH-1 is based on the piecewise- 
parabolic method (PPM) as described by Colella and 
Woodward The PPM scheme is a higher-order ex- 
tension of Godunov's method It uses parabolas as 
interpolation functions within each zone, and character- 
istics to determine the domains of dependence for zone 
interfaces. The average values of density, pressure and 
velocity within these domains are used as inputs to the 
Riemann problem between adjacent zones. The solution 
of the Riemann problem determines time averaged values 
of pressure and velocity, which in turn are used to com- 
pute hydrodynamical fluxes. In VH-l, each time step is 
reduced to three one-dimensional evolutionary "sweeps" 
via operator splitting. For each one-dimensional sweep, 
the fluid variables are evolved in Lagrangian coordinates 
and then remapped onto the original Eulerian grid. The 
order of the sweeps is cycled through all possible permu- 
tations. 

For the simulation presented here, I used a "Courant 
number" of 0.3. That is, the timestep was set to 0.3 times 
the maximum allowed by the CFL condition. This is a 
relatively small timestep for a PPM code. It was chosen 
to help minimize the artificial loss of angular momentum, 
discussed below. 

The Poisson equation for the gravitational potential is 
solved using multigrid methods [p^ . The finest grid has 
size 129'^, with grid points that lie at the corners of the 
128^ hydrodynamical zones. The boundary conditions 
are computed from a multipole expansion of the poten- 
tial that includes monopole, dipole, quadrupole, and oc- 
tupole terms. 

The PPM method is designed to evolve the fiuid mass 
density p, linear momentum density, and total energy 
density, and in the process to conserve the mass, lin- 
ear momentum, and total energy of the system. Unfor- 
tunately, for such a rapidly rotating star, the total en- 
ergy density E is dominated by kinetic energy density 
K. This is a problem because the internal energy den- 
sity pe = E ~ K is needed for the calculation of pressure 
from the equation of state (6). Since pe is a difference 
of large numbers, it is subject to large numerical errors. 
As long as the fluid flow is smooth, one can solve this 
problem by modifying the code so that internal energy, 
rather than total energy, is evolved. However, this mod- 
iflcation of the PPM algorithm leads to the wrong jump 
conditions at shock fronts. Thus, for the present simu- 
lation, a compromise was struck: in regions of smooth 
flow the code evolves internal energy, and in the neigh- 
borhood of shocks the code evolves total energy. Similar 
schemes have been described in Refs. [ |T7| , p^ . 



VH-1 uses the shock flattening algorithm described by 
Colella and Woodward (l4j as a dissipation mechanism. 
The algorithm determines a "flattening coefficient" / for 
each computational zone, which ranges from for smooth 
flow to a maximum of 0.5 in the presence of a strong 
shock. The flattening coefficient is also used to determine 
whether internal energy or total energy is used to update 
the fluid variables in each zone. Thus, internal energy 
is used when / is less than some threshold value, ft, 
while total energy is used for f > ft- In practice, the 
performance of the code was found to be insensitive to 
the value of the threshold ft, for values ranging from 0.1 
to 0.5. 

The ability of this "hybrid" code to switch from inter- 
nal energy to total energy at shock fronts was demon- 
strated on a variety of one-dimensional test problems. 
These include the standard Sod shock tube problem with 
density and pressure ratios of up to 1.0 x 10^, and strong 
standing shocks. Figures 2 and 3 show the results of the 
most extreme test, a standing shock wave with an up- 
stream Mach number of 1.0 x 10®. Results are shown for 
the hybrid code, a total energy code (which evolves total 
energy exclusively) and an internal energy code (which 
evolves internal energy exclusively) . The simulations use 
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FIG. 2. Logarithm of pressure for the strong standing 
shock test. The curves were obtained from the hybrid code 
(solid), total energy code (dashed) and internal energy code 
(dotted). 

100 zones to cover the computational domain, < a; < 1, 
and a Courant number of 0.6. The upstram fluid, in the 
region 0.5 < a; < 1, has density p = 1.0, pressure P = 1.0, 
and velocity u — —1.3 x 10*. The downstream fluid, in 
the region < x < 0.5, has density p = 4, pressure 
P = 1.25 X 10^6, and velocity u = -3.2 x lO'^. Ide- 
ally the discontinuities should be maintained a,t x — 0.5. 
Figures 2 and 3 show the pressure and density after 250 
timesteps. The CFL condition for this simulation is de- 
termined by the upstream velocity, so the upstream fluid 
travels a distance of ~ 0.006 in each timestep. 

Figure 2 shows that the total energy code produces 
extremely large errors in the upstream pressure. As ex- 
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FIG. 3. Density for the strong standing shock test. The 
curves are obtained from the three codes, as in Fig. 2. 

pected, Figs. 2 and 3 show that the results from the in- 
ternal energy code are quite poor. Altogether the hybrid 
code is the most successful at maintaining the proper 
pressure and density profiles. For this extreme test, the 
largest errors from the hybrid code appear in the down- 
stream fluid where the density and pressure differ from 
their ideal values by at most 11% and 1.6%, respectively. 
The hybrid code also lost about 2.7% of the total energy 
during the 250-timestep simulation. 

Because the present code conserves linear momentum, 
there is no erroneous center of mass motion due to nu- 
merical errors. On the other hand, the code does not 
conserve angular momentum. This is a fairly serious 
problem because the angular motion of the star is the 
central effect we wish to study. The amount of angu- 
lar momentum that is lost, or gained, through numerical 
error is found by computing the total angular momen- 
tum on the grid, Jgrid, and correcting for the amount of 
angular momentum that is lost (or gained) through the 
grid boundaries, Jiost- The angular momentum passing 
through the grid boundaries is obtained by computing 
the angular momentum that moves off the Eulerian grid 
during each Lagrangian time step. Ideally, the sum of the 
angular momentum on the grid and the angular momen- 
tum lost through the boundaries, Jgrid + Jiost, should be 
constant. As seen in the next section, over the course of 
the simulation (8000 time steps) about 8% of the star's 
angular momentum is lost through the boundaries. Dur- 
ing this same time, the total angular momentum drops 
by nearly 25%. Thus, numerical errors account for a drop 
of about 17% in the angular momentum of the star. 

In the Newtonian/quadrupole approximation, the 
gravitational wave signal hij is formed from linear com- 
binations of components of the second time derivative of 
the reduced quadrupole moment, lij. The relationship is 
hij = {2G/Rc'^)PIj^Im, where P^j^ is the projection oper- 
ator for transverse directions and R is the distance from 
the source to the observation point . Using the hydro- 
dynamical equations, the time derivatives that appear in 



i'ij can be eliminated |2^ so that, to within boundary 
terms, we have 

4 = STF J d^x2p[v^vj - x,dj'^>] . (7) 

Here, "STF" stands for the symmetric, trace free part of 
the expression that follows, Vi is the fluid velocity, dj — 
d/dxj is a spatial derivative, and <i> is the gravitational 
potential. 

The gravitational-wave signal is computed numerically 
by using expression (7) for lij. This calculation is subject 
to errors due to the finite size of the computational grid. 
Recall that during the course of the simulation, matter 
is expelled from the star. Some of this matter, a total of 
0.048 Mq, passes through the grid boundaries. As a con- 
sequence, the relationship hij = {2G / Rc^)PijIki is not 
correct, even in the Newtonian/quadrupole approxima- 
tion, because the lost matter is not included in the calcu- 
lation of lij . A second source of error stems from the fact 
that boundary terms were discarded in the derivation of 
Eq. (7). These boundary terms would vanish if the den- 
sity were always zero on the grid boundary. Fortunately, 
the matter that reaches the boundary has relatively low 
density and the errors that arise from the missing bound- 
ary terms in Eq. (7) are small. This has been verified by 
comparing the results from Eq. (7) with the results ob- 
tained by computing numerically the second time deriva- 
tive of the reduced quadrupole moment. The numerical 
derivatives were obtained by performing a least-squares 
fit to a quadratic using 5 or 7 consecutive data points. 
The second derivative at the central time is found from 
the curvature of the quadratic. This calculation is sub- 
ject to a certain amount of numerical noise, but otherwise 
the results agree quite closely with those obtained from 
Eq. (7). Since the matter that passes through the grid 
boundaries does not greatly effect the calculation of lij 
with Eq. (7), it is plausible to assume that the errors 
obtained by equating hij and {2G/Rd*')P}flke are also 
small. 

The frequency spectrum for the gravitational-wave 
signal is computed with the Lomb normalized peri- 
odogram [^. This method was chosen for its ability 
to handle the unevenly sampled data. In the Newto- 
nian/quadrupole approximation, the gravitational-wave 
luminosity is given by L — Glijl^^ / {5c^), and the ra- 
diated energy AE is the time integral of L |jl^. The 
calculation of the time derivative lij = dlij/dt must be 
handled with care since, otherwise, even small errors in 
L will accumulate in the integration over time. For this 
calculation, I use a Savitzky-Golay approach ||l^ . Specif- 
ically, the time derivative of the function lij in each time 
interval is computed from the slope of a quadratic poly- 
nomial that is obtained from a least-squares fit to several 
nearest neighbor data points. The number of points used 
is typically between 10 and 20. 

The evolution code was run on a Cray T90 vector 
computer, while the initial data, gravitational wave spec- 
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tnim, and gravitational wave Inminosity were computed 
on local workstations. The evolution required 8000 
timesteps and about 45 CPU hours. 



IV. DESCRIPTION OF THE EVOLUTION 

The star retains its predominantly circular, disk like 
shape for about 10 ms, which amounts to approximately 
6 rotation periods for the fluid near the center of the 
star. Over the next few milliseconds the bar mode grows 
rapidly, so that by 13 ms the star is highly elongated as 
shown in Fig 4. After a few more milliseconds the cen- 
tral regions of the star recircularize, and by 16 ms the 
inner density contours have nearly lost their bar shape 
as seen in Fig. 5. The cycle of bar formation and re- 
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FIG. 4. Contour plot of the density in the equatorial 
(X-Y) plane at 13.58 ms. The contour levels are 10^°g/cm^, 
10^^ g/c^l^ and lO" g/cm^. 

circularization repeats; by 19 ms the star again shows a 
strong bar-like shape, as shown in Fig. 6. The episodes of 
bar formation and recircularization continue with dimin- 
ishing amplitude and a period of about 6.5 ms. During 
each episode, the bar undergoes a little more than 1^ 
revolutions. Near the end of the simulation, at a time 
of ^ 50 ms, the oscillations become quite weak and the 
star settles into a bar-shaped configuration of modest 
strength. See Fig. 7. 

During the initial episode of rapid bar formation, mat- 
ter is thrown outward from the ends of the bar. The 
ejected matter forms long spiral arms, clearly visible in 
Figs. 4 and 5. Some of this matter is lost off the edge 
of the computational grid — between 14 and 20 ms the 
mass loss totals 0.040 Mq. The subsequent episodes of 



FIG. 5. Density contours as in Fig. 4, at 16.58 ms. 
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FIG. 6. Density contours as in Fig. 4, at 19.48 ms. 

bar formation are less violent, with only small amounts of 
ejected mass reaching the grid boundaries. The mass loss 
between 20 ms and the end of the simulation amounts to 
0.008 M©. 

When the star assumes its bar shape, the individual 
particles in the inner regions of the star circulate along 
trajectories that roughly coincide with the constant den- 
sity contours. For example, at the time shown in Fig. 4, 
13.58 ms, the bar is rotating about the center of mass 
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FIG. 7. Density contours as in Fig. 4, at the end of the 
simulation (52.57 ms). 



with a period of approximately 4 ms. Individual parti- 
cles flow along the lO^^g/cm^ density contour with fairly 
uniform speeds ranging from about 6.5 x lO^cm/s to 
about 7.5 X lO^cm/s. The density contour at lO^^g/cm^ 
has a perimeter length of approximately 1.7 x 10^ cm. 
Thus, the particles orbit along this contour with period 
~ 2.5 ms as the contour precesses about the star's center 
with period ~ 4 ms. 

The shape of the star can be quantified by Fourier 
analyzing the density in a circle in the equatorial plane 
1^,^. The Fourier coefhcients are defined by 
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dip p{ip)e'""f 



(8) 



where p{ip) is the density in the equatorial plane at an ar- 
bitrarily chosen distance of 2.0 x 10^ cm from the center of 
the star. Here, ip is the azimuthal angle. The amplitude 
of the m*-^ Fourier mode is defined by Cm — \J -^n + ^"im 
and the phase angle is 0^ = arctan(_B,„/74„j). Figure 8 
shows the natural logarithm of the ratio of the amplitude 
Cra and the average density Co for the m = 2 bar mode 
and the m = 4 mode. The m = 3 mode is small and is 
not displayed in Fig. 8. The coefficient C3/C0 remains 
less than 0.01 through most of the evolution, reaching 
peak values of ^ 0.02 in the late stages. 

The periodic growth and decay of the bar shape is 
clearly seen in the peaks and valleys of the graph in 
Fig. 8. The amplitudes grow exponentially for sev- 
eral milliseconds during the initial episode of bar for- 
mation. For the m — 2 bar mode, the growth rate near 
10ms is d\nC2/dt w 820/s. For the m = 4 mode, the 
growth rate is d\nC2/dt « 1400/s. In polytropic units 
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FIG. 8. Natural logarithm of the ratios C2/C0 and C4/C0 
of Fourier amplitudes. 



(r = 1.88 X 10~^s), the growth rates for the m = 2 
and TO = 4 modes are 0.015 and 0.026, respectively. For 
m — 2, this value differs by a few percent from the result 
0.0145 reported in Ref. pj. In units of the "dynamical 
time" to — [i?pq/(G'M)]^'*, where Rcq is the initial equa- 
torial radius, the growth rates are 0.50/ic and 0.86/t^) 
for the m = 2 and to = 4 modes. Comparing these re- 
sults to the values 0.55/iD and l.l/ic obtained from the 
highest resolution run in Ref. we find a fairly large 
discrepancy in the to = 4 case. This difference might be 
caused by numerical errors in the present code associated 
with the cartesian grid. 

Near the end of the simulation, the oscillations in the 
Fourier amplitudes have subsided and the star settles into 
a bar shape with strength C2/C0 « 0.18. The downward 
drift in C2/C0, visible in Fig. 8, is possibly due to nu- 
merical error as discussed below. 

The phase angles for the to = 2 and m — A modes 
are shown in Fig. 9. The eigenfrequencies are d(j)2/dt — 
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FIG. 9. Phase of the m = 2 and m = 4 Fourier modes 
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3.24/ms and d(h 4 / dt = 6.52/ms. The pattern speeds 
m~^d(t>m/dt are 1.62/nis for the m = 2 mode and 
1.63/ms for the m = 4 mode. These speeds, which equal 
0.99/ii3 and l.O/iu in dynamical time units, agree with 
the results obtained in Refs. |^J^ to within 1%. The ap- 
proximate equality of the pattern speeds shows that the 
m = 4 mode is a harmonic of the 171 — 2 bar mode. 

The stability parameter (3 undergoes fluctuations, 
dropping to a local minimum when the bar amplitude 
(^2 /Co is a local maximum, and rising to a local maxi- 
mum when the bar amplitude is a local minimum. This 
behavior is shown in Fig. 10. As the fluctuations dimin- 



rotation law of the type considered here, the dynamical 
m — 2 bar mode will not spontaneously grow unless (3 
exceeds a critical value near 0.27; however, once a bar 
shape is established, it can persist for many rotation pe- 
riods with (3 = 0.24 or less. 

As discussed in the previous section, the numerical 
code does not conserve angular momentum. In Fig. 11, 
the total angular momentum on the grid, Jgrid, is plotted 
as a function of time along with the difference between 
the initial angular momentum Jq — 6.98 x lO^^gcm^/s 
and the angular momentum Jiost that flows off the edge of 
the numerical grid. Ideally, these two curves should coin- 
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FIG. 10. The stability parameter /3 as a function of time. 

ish, (3 settles to a value of ~ 0.24 with a slight downward 
drift. Note that this value of (3 is well below the critical 
value 0.27 for the growth of the bar mode in a constant 
density star. 

The fact that the bar shape persists with /3 below 0.27 
might indicate that the threshold for bar formation in a 
7 — 5/3 star with rotation law given by Eq. (lb) is actu- 
ally 0.24 or less. This possibility was tested by evolving 
a second model star, obtained by modifying the freely 
specifiable data in Eq. (1) so that the central value of 
angular velocity is 3000/s rather than 4000/s. The pa- 
rameters describing this model are 



M ^ 


1.34 Mq , 


(9a) 


Rcq 


5.88 X 10^ cm , 


(9b) 


Rp = 


1.40 X 10^ cm , 


(9c) 


J = 


2.90 X lO'^^gcmVs , 


(9d) 


/3 = 


0.253 , 


(9e) 



so, in particular, the stability parameter is between 0.24 
and 0.27. This model was evolved for approximately 
23 ms, which equals about 11 rotation periods for the 
fluid near the center of the star. During this time the 
m = 2 bar mode showed no signs of growth. We are 
thus led to conclude that for a 7 = 5/3 fluid star with 
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FIG. 11. The angular momentum on the grid, Jgrid, and 
the difference Jo — Jioat between the initial angular momentum 
and the angular momentum lost at the boundaries of the grid. 

cide; the fact that Jgrid falls below Jq — Jiost indicates that 
numerical errors artificially remove angular momentum 
from the system. The results for Jiost were checked by 
computing the angular momentum flux from each of the 
saved data sets, which were generated every 50 timesteps. 
The numerically computed time derivative of Jiost is in 
very close agreement with the angular momentum flux. 

I have not yet determined the source of the numerical 
errors in angular momentum. Note, however, that the 
errors are not significant until the star begins to develop a 
noticeable bar shape, at about 10 ms. This suggests that 
the code has difficulty tracking the leading and trailing 
edges of the bar-shaped star, with their sharp density 
gradients, as the bar rotates in the x-y plane. 

Since numerical errors cause angular momentum to be 
lost, it is plausible to expect that rotational kinetic en- 
ergy is artificially lost as well. This might account for the 
downward drift in the stability parameter at late times, 
as seen in Fig. 10. In turn, one might expect a drop in 
the stability parameter to cause the star's bar shape to 
decay. Figure 8 shows that the Fourier amplitude C2 is 
fairly robust at late times with a slight downward trend 
on average. If angular momentum were not artificially 
lost, the simulation might show that after the fluctua- 
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tions cease the bar shape remains unchanged on a dy- 
namical time scale. 



V. GRAVITATIONAL WAVES 

In the Newtonian/quadrupole approximation, the 
gravitational wave amplitude for the plus polarization 
state, as measured at a distance R in the equatorial 
plane of the source, is 

c'^R 

-^h+{eq) = (cos^ ip - 2)1.^^ + (sin^ ip ~ 2)Iyy 

+2 cos ip sin ip I j;y . (10) 

Here, lij is the second time derivative of the reduced 
quadrupole moment and (p is the azimuthal angle of the 
observation point relative to the (arbitrary) x-axis of the 
source. Figure 12 shows the results for /i+ with R ~ 
20Mpc and (p ~ tt/2. (Recent estimates indicate that 



the XX and yy components satisiy Ixx + lyy ^ ^- These 
observations imply I^x ~ ~Iyy As a consequence, we see 
from Eqs. (10) and (11) that /i+(p) « 2/i+(eq). This is 
confirmed by the results of the simulation. As expected, 
the amplitude ft-x(p) is almost identical to /i+(p), just 
phase shifted by 45°. 

Although the reliability of the results for the 
gravitational- wave signal is rather poor, due both to the 
crudeness of the physical model and to numerical errors, 
it is still of some interest to consider the detectability of 
these waves. Thus, consider the characteristic amplitude 
of the source, approximated by 

Here, /S.E is the total energy radiated in gravitational 
waves and /c is the characteristic frequency. The 
gravitational-wave frequency is sharply peaked about 
490 Hz, as shown in Fig. 13. The total energy radi- 
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FIG. 12. The gravitational-wave amplitude fe+(eq) mea- 
sured at an observation point in the equatorial plane of the 
source at a distance of 20Mpc. 

the distance to the Virgo cluster is ~ 20.9 Mpc [Q.) In 
the equatorial plane, the cross polarization amplitude, 
^x(eq), is a linear combination of the components Ixz 
and lyz ■ Since the star is approximately symmetric under 
reflections in the equatorial plane, ft-x (eq) as measured in 
the equatorial plane is nearly zero. 

For an observation point above the north pole of the 
star, the plus polarization amplitude is 

c^R 

-— ft.+ (p) = (2 cos^ ip - l)Ixx + (2 sin^ -p - l)Iyy 

~\-Acospsmp Ixy ■ (11) 

Since the star is highly flattened and has little motion 
in the z-direction, the zz component of the second time 
derivative of the unreduced quadrupole moment, I^^, is 
very small. Furthermore, because the shape of the bar 
does not evolve rapidly compared to its rotation rate, 
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FIG. 13. The frequency spectrum of gravitational waves, 
computed with the Lomb normalized periodogram. 

ated in gravitational waves, as computed in the Newto- 
nian/quadrupole approximation up to the time the sim- 
ulation was halted, is AE « 5 x 1O~'*M0C^. From these 
data we find the characteristic amplitude of the source to 
be /ic ~ 4 X 10~^^ at a distance of 20 Mpc. This value for 
he is a lower bound, since the radiated energy AE will 
continue to increase as long as the bar persists. Also note 
that the result he w h^/n « 4 x 10^^^ coincides roughly 
with n w 25 wave cycles with amplitude h w 0.08 x 10~^^. 
This is consistent with the signal displayed in Fig. 12. 

The signal /ic ~ 4 x 10^^^ is probably too weak 
to be seen by the first generation of interferometric 
gravitational-wave detectors, which should have a sensi- 
tivity of about 10-20 at 500 Hz |6|j2|. The signal might 
be detectable by the planned advanced detectors, which 
are expected to have a sensitivity below 10"^^ at 500 Hz 
|]6|,p3t . For a source at a distance of ~ 25 kpc, comparable 
to the diameter of the Milky Way galaxy, the character- 
istic amplitude is /ic ~ 3 x 10~^^. This is a relatively 
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large signal that should be detected easily by the first 
generation of interferometers. 

If the simulation is scaled as described in section II, 
the amplitude as measured at a fixed (unsealed) distance 
R changes by a factor of r^/(/iA^) and the frequency 
changes by a factor of t. For example, with A = 1/2, fi = 
2, and r — 1/4, the amplitude h+{eq) at R — 20Mpc is 
reduced from that shown in Fig. 12 by a factor of 8. The 
characteristic amplitude he is also reduced by a factor of 
8, to about 5 x 10^^^ at 20 Mpc, and the peak frequency 
occurs at ^ 120 Hz. This signal is somewhat below the 
expected sensitivity of the advanced gravitational-wave 
detectors, which is approximately 10~^^ at 120 Hz P,p3[. 
At a distance of 25 kpc, the characteristic amplitude is 
almost 4 x 10~^°, well above the expected sensitivity of 
~ 4 X 10~^^ for the first generation detectors. 



the first generation interferometers. 
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VI. CONCLUSIONS 



The simulation presented here shows that a rapidly 
rotating 7 = 5/3 fluid star with rotation law Eq. (lb) is 
dynamically unstable to bar formation, and that the bar 
shape, once established, is relatively long lived. Over the 
final 10 or 20 ms of the simulation, the star maintains 
a modest but fairly persistent bar shape with stability 
parameter f3 equal to about 0.24. It is possible that, 
for 7 = 5/3 stars with certain initial angular velocity 
profiles, the threshold for bar formation is actually 0.24 
or less. However, further numerical study showed that 
an axisymmetric star with (3 = 0.25 initially is stable 
against growth of the bar mode. Thus, for the models 
considered here, the dynamical bar mode requires a value 
of P near 0.27 to grow but the bar shape can persist 
for many bar-rotation periods with /3 = 0.24 or less. 
Imamura, Durisen, and Pickett suggest that as long as /3 
is greater than the critical value for growth of the secular 
instability, about 0.14, a bar-shaped star is dynamically 
stable against forming an axisymmetric disk The 
results here are consistent with this conjecture. Recent 
numerical studies by New, Centrella, and Tohline |^ also 
indicate that the bar shape is long lived. 

The strength of the gravitational-wave signal emitted 
by an initially axisymmetric, dynamically unstable star 
depends on the rotating star's length and time scales. 
With the numerical model scaled to represent a low- 
density neutron star with soft equation of state, the char- 
acteristic amplitude is /ic « 4 x 10~^^ at a distance of 
20 Mpc. With the model scaled to represent a stellar 
core that has partially collapsed and is prevented from 
contracting to nuclear densities by centrifugal forces, the 
characteristic amplitude is /ic ~ 5 x lO^^'^ at 20 Mpc. 
These signals are probably too weak to be detected by 
the first generation of interferometers, although the neu- 
tron star signal might be detectable by the planned ad- 
vanced interferometers. At a distance of 25 kpc, both 
signals are quite large and should be detected easily by 
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APPENDIX 



The technique used in this work to solve the equations of hydrostatic equihbrium is similar to the self-consistent 
field method, first developed by Ostriker and Mark and later refined by Hachisu (25|. The self-consistent field 
method is based on the integral form of the Euler equation which reads, for the equation of state P = Kp^ [p4[[25|], 



p^-^ + * + $ = C . (13) 



7-1 



In this equation, ^E" — — Jdrruj^ is the "rotational potential" derived from the angular velocity <i> is the 

gravitational potential, and C is an integration constant. In the self-consistent field method, one begins with an 
initial guess for the density p, solves the Poisson equation for $, then uses Eq. (13) to obtain a corrected density 
distribution. The process is then iterated. Typically, the Poisson equation is solved by expanding the potential 
in terms of Legendre polynomials, and the constant C appearing in Eq. (13) is determined by specifying certain 
properties at the star's surface. For example, Hachisu fixes the "axis ratio" Rp/Rcq [p5| . 

For the method used in this paper, the constant C is determined by properties at the center of the star. The 
rotational potential is chosen to vanish at the star's center, ^'(0) — 0, and the constant C is written as 

C - jKp2-^ + $c . (14) 

Here, the subscript c denotes the center of the star. The central density pc is specified as input data, along with the 
rotation law Lu{r). The initial guess for p is chosen to be a static, spherically symmetric polytrope with central density 
Pc, obtained by solving numerically the Lane-Emden equation. The potential and its central value, are obtained 
from the Poisson equation using a multigrid algorithm, as discussed below. Note that with this scheme the center of 
the star must have nonzero density. This precludes the possibility of generating stars with toroidal surfaces. 
Equation (13) can be written as F = where the function F is defined by 

F = J^p^-^ + + q>-C (15) 

and the constant C is given by Eq. (14). The corrections to the density are obtained from a Newton-Raphson 
algorithm applied to F = 0. Note that F{x) is a nonlocal function of p{x), where x labels points in space. The 
nonlocal dependence of F{x) on p{x) enters through the gravitational potential $(a;). One might expect that the 
value of <I> at the point x is insensitive to the density p{x) at the same point x, since the potential $ is determined 
primarily by the global properties of the star, such as its mass. With this observation in mind, we compute the 
variation of F at point x with respect to p at point x by dropping the nonlocal terms. This leads to the approximate 
result 6F « p'^~^5p. Then setting F + 5F = 0, we find the density correction 

5p^^p^~^{C^^~^)--^p , (16) 

7 ~ 1 

to be used for each Newton-Raphson iteration. 

I have not constructed a mathematical argument to justify the above assumption, namely, that the nonlocal terms 
in F can be ignored in computing 5F . However, this same assumption is implicit in the usual self-consistent field 
method. As described in Refs. p^ , p5[ , for the self-consistent field method the density is updated by solving Eq. (13) 
for p, ignoring the p-dependence in $. That is, the change in p is 



'"-P, (17) 

where $ is the solution of the Poisson equation. Formula (17) can be justified by the same reasoning that led to 
Eq. (16), but with the following difference: enthalpy H = p'^~^ / {'^ — 1), rather than density p, is treated as the 
independent variable. The variation of F with respect to ff, ignoring the $ dependence on H , is 5F w 5H . Setting 
F + 5F = leads to a new enthalpy distribution, -ffnow ^ H ~ F. Solving this equation for the new density Pnowj we 
find the density correction 5p = pncw — p of Eq. (17). 

The initial data code solves the equations of hydrostatic equilibrium at the zone centers of a 513 x 513 grid in the 
r-z plane of a cylindrical coordinate system. The symmetry axis concides with the edge of the first column of zones. 
The density p is computed on the inner 512 x 511 zones. One layer of zones at the top, bottom, and side opposite the 
symmetry axis are reserved for setting boundary conditions for the gravitational potential The boundary values 



5p- 
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(C-$-*) 
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are obtained by expanding the potential in a multipole series that includes monopole and quadrupole terms. (The 
dipole and octupole moments vanish due to rotational symmetry and reflection symmery about the equatorial plane.) 
The Poisson equation is solved with a multigrid algorithm, adapted from the code described in Ref. ||l^ to handle the 
boundary conditions and the cylindrical geometry. 

The initial data described in Sec. II was generated by a code that relies on the Newton-Raphson iteration scheme 
(16) and the multigrid Poisson solver described above. The results of three tests of this code (JDB) are shown in 
Table I. For comparison, the results obtained by Hachisu ||2^ are also listed. The values reported in the table are 
scaled as in Hachisu's self-consistent field method, with Newton's constant G, the equatorial radius Rcq, and the 
maximum density Pmax equal to unity. Hats denote these scaled quantities. Also, 11 denotes the volume integral of 
pressure, / PdV. 





test 


Rp 


M 


J 


f 


-W 


3n 


p 

^ max 


V 




JDB 


r 


0.6667 


0.3288 


0.02575 


0.006641 


0.1164 


0.1031 


0.2044 


1.8 X 10- 


4 


Hachisu 


r 


0.667 


0.328 


0.0257 


0.00663 


0.116 


0.103 


0.204 






JDB 


V 


0.3332 


0.6413 


0.1378 


0.06392 


0.3733 


0.2454 


0.2020 


4.7 X 10^ 


5 


Hachisu 


V 


0.333 


0.639 


0.137 


0.0638 


0.372 


0.244 


0.202 






JDB 


3 


0.1662 


0.8419 


0.1036 


0.04559 


0.5982 


0.5070 


0.3272 


1.6 X 10" 


5 


Hachisu 


j 


0.167 


0.843 


0.104 


0.0456 


0.599 


0.508 


0.327 







The "r" test assumes a rigid rotation law, Lu{r) — ujq, with the scaled angular velocity given by luq — 0.266. The 
"v" test assumes the "w-constant" rotation law, uj{r) — va/\uP + r^, with Vq — 0.215 and d = 0.100. The "j" test 
assumes the "j-constant" rotation law, uj{r) = jo/{d^ + r^), with jg = 0.0176 and d = 0.100. For each test, 7 = 5/3. 
For Hachisu's data the virial theorem diagnostic V is typically less than a few times 10"'* ||2^. The test results show 
that the present code and Hachisu's code differ by less than 1% in every case. 

Subsequent testing has shown that the density update in Eq. (16) is no better, in terms of accuracy and convergence 
rate, than the usual self-consistent field method update (17). After 30 Newton-Raphson iterations, the results 
obtained from the two update formulas agree to a few parts in 10~^, at worst. For both update formulas, the errors 
(as compared to the 30-iteration results) were at most a few tenths of a percent after 10 iterations. 
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